1. Wczytanie zestawu danych punktowych oraz nadanie im ukladu wspolrzednych ETRS 1989 Poland CS2000 Zone 7 - EPSG:2178, nastepnie zapisanie ich do nowego pliku csv - zestaw1_epsg_2178.csv.
library(sp)
data_points <- read.csv('zestaw1.csv', colClasses = c('numeric', 'numeric'))
head(data_points)
## Long Lat
## 1 19.89514 50.07248
## 2 20.03497 50.07282
## 3 19.93090 50.04349
## 4 19.86676 50.06638
## 5 19.93647 50.05540
## 6 19.93338 50.06552
dim(data_points)
## [1] 2000 2
coord <- SpatialPoints(cbind(data_points$Long, data_points$Lat),
proj4string = CRS('+proj=longlat'))
class(coord)
## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"
length(coord)
## [1] 2000
head(coord)
## SpatialPoints:
## coords.x1 coords.x2
## [1,] 19.89514 50.07248
## [2,] 20.03497 50.07282
## [3,] 19.93090 50.04349
## [4,] 19.86676 50.06638
## [5,] 19.93647 50.05540
## [6,] 19.93338 50.06552
## Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84
## +no_defs
coordUTM <- spTransform(coord, CRS('+init=epsg:2178'))
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum European_Terrestrial_Reference_System_1989 in
## Proj4 definition
head(coordUTM)
## SpatialPoints:
## coords.x1 coords.x2
## [1,] 7420912 5549067
## [2,] 7430922 5548966
## [3,] 7423426 5545805
## [4,] 7418870 5548418
## [5,] 7423844 5547124
## [6,] 7423639 5548253
## Coordinate Reference System (CRS) arguments: +proj=tmerc +lat_0=0
## +lon_0=21 +k=0.999923 +x_0=7500000 +y_0=0 +ellps=GRS80 +units=m
## +no_defs
length(coordUTM)
## [1] 2000
#uklad CRS pomyslnie zmieniony z WGS84 na EPSG:2178
coordUTM_df <- as.data.frame(coordUTM) #col.names nie dziala - zbadac, z jakiego powodu
colnames(coordUTM_df) <- c('Lon', 'Lat')
head(coordUTM_df)
## Lon Lat
## 1 7420912 5549067
## 2 7430922 5548966
## 3 7423426 5545805
## 4 7418870 5548418
## 5 7423844 5547124
## 6 7423639 5548253
dim(coordUTM_df) #dlugosc zgadza sie z dlugoscia wczytanych danych
## [1] 2000 2
write.csv(coordUTM_df, "zestaw1_epsg_2178.csv", row.names = FALSE)
2. Wczytanie danych w ukladzie wspolrzednych EPSG:2178 oraz pliku shapefile zawierajacego osiedla.
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.5 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.4 ✓ stringr 1.4.0
## ✓ readr 2.0.2 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(ggplot2)
library(raster)
##
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
##
## select
data <- read.csv('zestaw1_epsg_2178.csv')
head(data)
## Lon Lat
## 1 7420912 5549067
## 2 7430922 5548966
## 3 7423426 5545805
## 4 7418870 5548418
## 5 7423844 5547124
## 6 7423639 5548253
dim(data)
## [1] 2000 2
str(data) #2 zmienne - Lon oraz Lat numerical
## 'data.frame': 2000 obs. of 2 variables:
## $ Lon: num 7420912 7430922 7423426 7418870 7423844 ...
## $ Lat: num 5549067 5548966 5545805 5548418 5547124 ...
sum(is.na(data)) #brak pustych danych
## [1] 0
qplot(data$Lon, geom='histogram')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
districts <- shapefile('osiedla.shp')
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS = dumpSRS, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=tmerc +lat_0=0 +lon_0=21 +k=0.999923 +x_0=7500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
## but +towgs84= values preserved
Wczytane dane maja rozmiar 2000 x 2, zaiweraja dlugosc i szerokosc geograficzna - typ numerical, zapisane sa w data.frame o nazwie data, brak pustych danych.
3. Prezentacja mapy dzielnic/osiedl w Krakowie z nalozonymi na nia wczytanymi punktami.
map_cracow <- ggplot() +
geom_polygon(data = districts, aes(x = long, y = lat, group = group), show.legend = FALSE, color = "white", fill = "darkgrey") + coord_fixed()
## Regions defined for each Polygons
summary(districts)
## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x 7413437 7443955
## y 5537344 5555031
## Is projected: TRUE
## proj4string :
## [+proj=tmerc +lat_0=0 +lon_0=21 +k=0.999923 +x_0=7500000 +y_0=0
## +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs]
## Data attributes:
## GESTOSC_ZA ID IL_K_MOBIL IL_K_NIEMO
## Min. : 0 Length:141 Length:141 Length:141
## 1st Qu.: 293 Class :character Class :character Class :character
## Median : 1625 Mode :character Mode :character Mode :character
## Mean : 5671
## 3rd Qu.: 7722
## Max. :30031
## IL_K_POPRO IL_K_PROD IL_K_PRZED IL_L_MOBIL
## Length:141 Length:141 Length:141 Length:141
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## IL_L_NIEMO IL_L_POPRO IL_L_PROD IL_L_PRZED
## Length:141 Length:141 Length:141 Length:141
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## IL_M_MOBIL IL_M_NIEMO IL_M_POPRO IL_M_PROD
## Length:141 Length:141 Length:141 Length:141
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## IL_M_PRZED ILOSC_K_10 ILOSC_KOBI ILOSC_MEZC
## Length:141 Min. : 0.0 Length:141 Length:141
## Class :character 1st Qu.:103.3 Class :character Class :character
## Mode :character Median :110.8 Mode :character Mode :character
## Mean :113.2
## 3rd Qu.:119.6
## Max. :488.7
## ILOSC_MIES KOD MAPID MSLINK
## Length:141 Length:141 Length:141 Length:141
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## NAZWA_JEDN NR_JEDN_UR OZN_JEDN_U POWIERZCHN
## Length:141 Length:141 Length:141 Min. : 83366
## Class :character Class :character Class :character 1st Qu.: 892431
## Mode :character Mode :character Mode :character Median : 1638151
## Mean : 2318933
## 3rd Qu.: 3018523
## Max. :11796876
## WSP_X WSP_Y
## Min. :21980 Min. :-301839
## 1st Qu.:28779 1st Qu.:-295035
## Median :31681 Median :-292382
## Mean :31075 Mean :-291027
## 3rd Qu.:33641 3rd Qu.:-286706
## Max. :37118 Max. :-275075
map_cracow + geom_point(data=data, aes(x=Lon, y=Lat), alpha=0.8, size = 1.2, colour="darkred") +ggtitle("Wykroczenia na mapie Krakowa\n z podzialem na dzielnice") +
theme(plot.title = element_text(hjust = 0.5)) +
labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")
4. Wydzielenie klastrow o zwiekszonej intensywnosci wystepowania wykroczen przy uzyciu algorytmow grupowania gestosciowego DBSCAN/ HBDSCAN / OPTICS.
4a
• DBSCAN (density-based spatial clustering of applications with noise) - nieliniowy algorytm uczenia sie bez nadzoru, ktory wykorzystuje technikę osiągalności i łączności gęstości. DBSCAN dzieli dane na grupy punktow majace wspolne cechy charakterystyki lub skupienia.
• Zasada działania algorytmu:
1) Losowo wybierany jest punkt p
2) Pobierane sa wszystkie punkty, ktorch gestosc jest osiagalna wzgledem punktu p na podstawue maksymalnego promienia sasiedztwa (EPS) oraz minimalnej liczby punktow w sasiedzwie EPS (MinPts)
3) Dla kluczowych puntkow p tworzony zostaje klaster, w innym przypadku punkt p zostaje zaklasyfikowany jako punkt odstajacy lub szum.
• Klasyfikacja punktów:
1) Centralny (core point) - przynajmniej MinPts liczba punktow (wlacznie z samym punktem p) znajduje sie w otoczeniu p z promieniem rownym EPS
2) Graniczny (border point) - jezeli punkt jest osiagalny z punktu centralnego i w otoczeniu znajduje sie liczba punktow mniejsza, niz MinPts
3) Odstajacy - jezeli punkt nie jest ani granicznym ani centralnym zotaje uznany za odstajacy
• Plusy DBSCAN:
- dobre dzialanie przy dowolnych ksztaltach klastrow
- odporny na wartosci odstajace, jest w stanie je wykryc
- nie trzeba wczesniej z gory okreslac liczby klastrow
• Minusy DBSCAN:
- charakterystyka tworzonych klastrow jest okreslana przez parametry MinPts oraz EPS przez co jezeli ciezko jest czasami utworzyc klastry o znaczacych roznicach w gestosci
- czasami wyznaczenie odpowiedniej odleglosci sasiedztwa EPS nie nalezy do latwych i wymaga dodatkowej wiedzy
• Zastosowanie DBSCAN z parametrami odpowiednio MinPts i EPS rownymi:
5/10, 10/50, 10/100, 10/200
library(dbscan)
dbscan_clusters1 <- dbscan(data, minPts = 5, eps = 10)
dbscan_clusters1
## DBSCAN clustering for 2000 objects.
## Parameters: eps = 10, minPts = 5
## The clustering contains 50 cluster(s) and 1569 noise points.
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
## 1569 7 8 10 13 5 7 11 5 8 10 11 6 14 9 12
## 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
## 11 11 14 13 18 5 20 7 9 7 5 10 8 7 5 9
## 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
## 6 6 5 8 6 6 6 13 5 16 7 11 9 5 6 5
## 48 49 50
## 5 6 5
##
## Available fields: cluster, eps, minPts
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = dbscan_clusters1$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("DBSCAN, minPts=5, eps=10") +
theme(plot.title = element_text(hjust = 0.5)) +
labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")
dbscan_clusters2 <- dbscan(data, minPts = 10, eps = 50)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = dbscan_clusters2$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("DBSCAN, minPts=10, eps=50") +
theme(plot.title = element_text(hjust = 0.5)) +
labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")
dbscan_clusters3 <- dbscan(data, minPts = 10, eps = 100)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = dbscan_clusters3$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("DBSCAN, minPts=10, eps=100") +
theme(plot.title = element_text(hjust = 0.5)) +
labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")
dbscan_clusters4 <- dbscan(data, minPts = 10, eps = 200)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color =dbscan_clusters4$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("DBSCAN, minPts=10, eps=200") +
theme(plot.title = element_text(hjust = 0.5)) +
labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")
• Wnioski z DBSCAN:
1) Dla minPts = 5 i EPS = 10 bardzo mala ilosc klastrow, najwieksze skupisko w rejonie dzielnicy Stare Miasto, nastepnie mozemy zauwazyc klastry o zdecydowanie mniejszej gestosci w dzielniach Pradnik Bialy, Pradnik Czerwony, Mistrzejowice oraz Bienczyce. Szczatkowe ilosci punktow wystepuja rowniez na poludniu w dzielnicy Swoszowice.
2) Dla minPts = 10 i EPS = 50 zaczynamy zauwazac wyrazne klastry w okolicy dzielnicy Starego Miasta oraz pojedyncze wystepujace na terenie Nowej Huty/Wzgorz Krzeslawickich.
3) Dla minPts = 10 i EPS = 100 coraz wiecej klastrow o duzej gestosci w dzielnicy Starego Miasta, Grzegorzki, Czyzyny oraz Nowej Huty.
4) Dla minPts = 10 i EPS = 200 najwieksze, bardzo wyrazne skupisko klastrow w dzielnicy Stare Miasto, nastepnie mozemy zauwazyc wieksze klastry punktow w dzielnicy Krowodrza, Pradnik Bialy, Pradnik Czerwony, Mistrzejowice, Grzegorzki. Mniejsze klastry znajdziemy w Czyzynach, Wzgorzach Krzeslawickich oraz Nowej Hucie. Pojedynczy klaster wystepuje w Podgorzu Duchackim.
4b
• HDBSCAN (hierarchical density-based spatial clustering) - bardziej zaawansowana wersja algorytmu DBSCAN. Algorytm wykorzystuje podejscie oparte na gestosci - zamiast szukac klastrow posiadajacych okreslony ksztalt szuka regionow danych, ktore sa gestsze od otaczajacej je przestrzeni.
• Plusy HDBSCAN:
- dla danych o mocno zroznicowanej gestosci jest lepszy oraz szybszy w porownaniu do DBSCAN
- HDBSCAN w czasie dzialania odrzuca male odrosty punktow zachowujac najwieksze klastry okreslone przez parametr minimalnego rozmiaru klastra - dzieki temu rowniez dendogram algorytmu HDBSCAN jest mniej skomplikowany
• Minusy HDBSCAN:
- mimo szybszego czasu dzialania jest duzo ciezszy w zrozumieniu, przez wiele roznych operacji, ktore algorytm wykonuje w czasie dzialania
• Zastosowanie HDBSCAN z parametrem minPts wynoszacym odpowiednio:
10, 50, 100, 200
hdbscan_cluster1 <- hdbscan(data, minPts = 10)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = hdbscan_cluster1$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("HDBSCAN, minPts=10") +
theme(plot.title = element_text(hjust = 0.5)) +
labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")
hdbscan_cluster2 <- hdbscan(data, minPts = 50)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = hdbscan_cluster2$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("HDBSCAN, minPts=50") +
theme(plot.title = element_text(hjust = 0.5)) +
labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")
hdbscan_cluster3 <- hdbscan(data, minPts = 100)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = hdbscan_cluster3$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("HDBSCAN, minPts=100") +
theme(plot.title = element_text(hjust = 0.5)) +
labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")
hdbscan_cluster4 <- hdbscan(data, minPts = 150)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = hdbscan_cluster4$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("HDBSCAN, minPts=150") +
theme(plot.title = element_text(hjust = 0.5)) +
labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")
• Wnioski z HDBSCAN:
1) Przy parametrze minPts = 10 dominuje duzy klaster obejmujacy dzielnice Stare Miasto, a takze poludniowe dzielnice sasiadujace takie jak Podgorze,
Lagiewniki, wschodnie Grzegorzki, Czyzyny, Bienczyce, Mistrzejowice, oraz polnocne Pradnik Bialy i Czerwony.
2) Przy parametrze minPts = 50 widac juz dwa klastry, pierwszy najwiekszy obejmujacy Stare Miasto i okolice, drugi w obszarze Mistrzejowic, Bienczyce oraz Nowa Hute.
3) Dla parametru minPts = 100 brak wiekszych zmian w porownaniu do wyniku z minPts = 50.
4) Dla parametru minPts = 150 brak klastrow.
optics_out1 <- optics(data, minPts = 10)
optics_cluster1 <- extractDBSCAN(optics_out1, eps_cl = 10)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = optics_cluster1$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("OPTICS, minPts=10, eps_cl=10") +
theme(plot.title = element_text(hjust = 0.5)) +
labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")
optics_out2 <- optics(data, minPts = 10)
optics_cluster2 <- extractDBSCAN(optics_out2, eps_cl = 50)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = optics_cluster2$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("OPTICS, minPts=10, eps_cl=50") +
theme(plot.title = element_text(hjust = 0.5)) +
labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")
optics_out3 <- optics(data, minPts = 10)
optics_cluster3 <- extractDBSCAN(optics_out3, eps_cl = 100)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = optics_cluster3$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("OPTICS, minPts=10, eps_cl=100") +
theme(plot.title = element_text(hjust = 0.5)) +
labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")
optics_out4 <- optics(data, minPts = 10)
optics_cluster4 <- extractDBSCAN(optics_out4, eps_cl = 650)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = optics_cluster4$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("OPTICS, minPts=10, eps_cl=650") +
theme(plot.title = element_text(hjust = 0.5)) +
labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")
• OPTICS (ordering points to identify the clustering structure) - algorytm wykorzystujacy drzewa kd, stworzony w celu wyeliminowania wad algorytmu DBSCAN. OPTICS w trakcie dzialania dynamicznie rozszerza promien wyszukiwania dookola kazdego przypadku, zamiast ustalania z gory okreslonej wartosci promienia. Po zakonczeniu dzialania, a wie po przejsciu wszystkich przypadkow zwracana jest kolejnosc przetwarzania - "odwiedzin" oraz wynik osiagalnosci.
• Zalety:
- brak narzucania z gory wielkosci promienia
- nie wymaga parametrow gestosci
- kolejnosc klastrowania moze byc uzyteczna przy wyodrebnianiu informacji o klastrowaniu
• Wady:
- algorytm nie radzi sobie z wielkowymiarowymi danymi
- tworzy tylko porzadek klastrowy